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Abstract 


2  A  variational  method  is  used  to  estimate  wave-affected  parameters  in  a 

3  two-equation  turbulence  model  with  assimilating  the  temperature  data  into  an  ocean 

4  boundary  layer  model.  Enhancement  of  turbulent  kinetic  energy  dissipation  due  to 

5  breaking  waves  is  considered.  The  Mellor-Yamada  2.5  turbulence  closure  scheme 

6  (MY-2.5)  with  the  two  uncertain  wave-affected  parameters  (wave  energy  factor  a  and 

7  Chamock  coefficient  p)  is  selected  as  the  two-equation  turbulence  model  for  this 

8  study.  Two  types  of  experiments  are  conducted.  First,  within  an  identical  synthetic 

9  experiment  framework,  the  upper  layer  temperature  “observations”  in  summer 

10  generated  by  a  “truth”  model  are  assimilated  into  a  biased  simulation  model  to 

11  investigate  if  (a,  p)  can  be  successfully  estimated  using  the  variational  method. 

12  Second,  real  temperature  profiles  from  the  Ocean  Weather  Station  Papa  are 

13  assimilated  into  the  biased  simulation  model  to  obtain  the  optimal  wave-affected 

14  parameters.  With  the  optimally-estimated  parameters,  the  upper  layer  temperature  can 

15  be  well  predicted.  Furthermore,  the  horizontal  distribution  of  the  wave-affected 

16  parameters  employed  in  a  high  order  turbulence  closure  scheme  can  be  estimated 

17  optimally  by  using  the  four-dimensional  variational  method  that  assimilates  the  upper 

18  layer  available  temperature  data  into  an  ocean  general  circulation  model. 
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1.  Introduction 


Observations  (Kitaigorodskii  et  al.  1983;  Thorpe  1984,  1992;  Anis  and  Mourn  1992; 
Terray  el  al.  1996;  Drennan  et  al.  1996)  show  that  the  dissipation  of  turbulent  kinetic 
energy  (TKE)  is  enhanced  greatly  near  the  sea  surface  due  to  increasing  shear  by 
surface  gravity  waves  under  non-breaking  (including  Langmuir  circulation)  and 
breaking  waves.  The  mixing  induced  by  non-breaking  waves  directly  affects  or 
influences  the  upper-ocean  mixing  down  to  depths  of  the  order  of  100  m.  With  a  wave 
amplitude-based  Reynolds  number  (Re),  an  empirically  determined  critical  value  (Recr) 
is  used  to  identify  if  the  turbulence  is  generated  by  waves  (Re  >  Recr)  or  not  (Re  < 
Recr);  and  to  determine  a  depth  of  the  upper  ocean  mixed  layer  from  Re  =  Recr. 
Decrease  of  Re  with  depth  confirms  that  within  that  depth  the  turbulence  is  generated 
by  orbital  movement  of  surface  gravity  waves;  and  below  that  depth  there  is  no 
wave-induced  turbulence  (. Babanin  2006;  Babanin  and  Haus  2009). 

The  breaking  wave  induced  mixing  has  been  broadly  implemented  into  ocean 
circulation  and  mixing  models  (e.g.,  Craig  1994).  On  the  base  of  the  observational 
evidences  on  the  surface  wave  breaking  (Osborn  et  al.  1992;  Agrawal  et  al.  1992), 
Terray  el  al.  (1996)  suggested  a  three-layer  structure:  The  first  layer  (from  the  surface) 
is  a  wave-enhanced  layer  with  the  depth  on  the  same  order  as  the  significant  wave 
height,  and  the  energy  dissipation  rate  proportional  to  z“3  ( z  denotes  the  vertical 
distance  from  the  sea  surface),  which  is  twice  faster  than  the  classical  wall  layer 
dissipation.  The  second  layer  is  the  transition  layer  below  the  breaking  zone  (depth 
about  6z0,  z0  the  surface  roughness  length)  (Craig  and  Banner  1994),  with  the 
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1  energy  dissipation  rate  proportional  to  z“2 .  The  third  layer  is  the  classic  wall  layer 

2  with  the  energy  dissipation  rate  proportional  to  depth  z"1. 

3  To  model  the  wave-breaking  enhanced  turbulence  near  the  sea  surface  layer, 

4  Craig  and  Banner  (1994,  1996)  imposed  a  surface  diffusion  boundary  condition  on  the 

5  turbulent  kinetic  energy  equation  (hereafter,  CB  boundary  condition)  in  the 

6  Mellor-Yamada  (MY)  turbulence  closure  model  (1982).  Burchard  (2001b)  simulated  a 

7  wave-enhanced  layer  under  breaking  surface  waves  with  a  two-equation  turbulence 

8  model  including  the  CB  boundary  condition.  Mellor  and  Blumberg  (2004)  developed 

9  a  wave-enhanced  parameterization  scheme  with  the  CB  boundary  condition  to 

10  overcome  a  weakness  of  the  MY  turbulence  closure  model  that  produces  a  shallower 

1 1  surface  boundary  layer  and  higher  surface  temperature  during  summertime  warming 

12  in  comparison  to  the  observations  (Martin,  1985).  Zhang  et  al  (2011,  2012)  identified 

13  the  effect  of  breaking  surface  waves  on  upper  ocean  boundary  layer  deepening  in  the 

14  Yellow  Sea  in  summer  utilizing  the  Princeton  ocean  model  (POMgcs,  Ezer  et  al., 

15  2004).  A  well-mixed  temperature  surface  layer  in  the  Yellow  Sea  can  be  reconstructed 

16  successfully  when  the  breaking  wave  enhanced  turbulent  mixing  are  considered. 

17  In  addition  to  the  wave  breaking,  other  wave-related  processes  are  also  important 

18  in  modulating  the  upper  mixed  layer,  such  as  the  non-breaking  wave  (Babanin  et  al, 

19  2009)  and  the  Langmuir  circulation  (Stephen  et  al,  2012).  Some  studies  indicate  that 

20  the  effect  of  wave  breaking  on  the  upper-level  turbulence  is  significant  within  the 

21  depth  comparable  to  the  wave  height  (Terray  et  al.,  1996;  Babanin  et  al,  2005). 

22  However,  for  a  deeper  mixed  layer  over  100  m  depth,  the  impact  of  wave  breaking 
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1  would  be  small  and  the  effect  of  Langmuir  circulation  and  non-breaking  wave 

2  becomes  important  (Babanin,  2005). 

3  Uncertain  wave-affected  parameters  exist  in  modeling  wave-induced  turbulence 

4  (non-breaking  or  breaking  waves)  such  as  critical  value  of  the  wave  Renolds  number 

5  (Recr)  in  non-breaking  waves  and  wave  energy  factor  (a)  and  Chamock  coefficient  (fi) 

6  in  breaking  waves.  These  parameters  are  usually  determined  empirically  or  adjusted 

7  artificially.  Studies  have  shown  successful  parameter  estimation  with  a  dynamical 

8  model  using  variational  optimal  control  techniques  (Derber,  1987;  Le  Dimet  and 

9  Talagrand,  1986).  For  example,  Yu  and  O'Brien  (1991)  used  the  variational  method  to 

10  assimilate  meteorological  and  oceanographic  observations  into  a  one-dimensional 

11  oceanic  Ekman  layer  model,  to  estimate  the  drag  coefficient  and  the  oceanic  eddy 

12  viscosity  profde,  and  to  investigate  the  effect  of  initial  condition  on  the  variational 

13  parameter  estimation.  Zhang  et  al  (2003)  showed  the  capability  of  4D-Variational 

14  method  (4D-VAR)  in  estimating  uncertain  parameters  in  numerical  models.  Peng  et  al 

15  (2006)  developed  a  tangent  linear  model  and  an  adjoint  model  of  three-dimensional 

16  POM  to  construct  a  4D-VAR  algorithm  for  coastal  ocean  prediction.  Effective  error 

17  correction  was  found  in  initial  conditions  and  wind  stress  in  the  storm  surge 

18  simulation  (Peng  et  al,  2007),  and  the  drag  coefficient  was  estimated  in  the  storm 

19  surge  prediction  using  the  adjoint  model  of  the  three-dimensional  POM  [Peng  et  al 

20  (2012)].  Peng  et  al  (2006)  also  pointed  out  that  it  is  still  an  open  issue  as  to  whether  it 

21  is  meaningful  to  linearize  the  turbulence  closure  scheme  in  an  atmospheric  or  oceanic 

22  model  due  to  the  high  nonlinearity  and  discontinuity  of  the  vertical  turbulence.  The 
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1  nonphysical  noise  might  be  produced,  and  thus  lead  to  numerical  instability  during  the 

2  process  of  linearizing  the  turbulence  closure  scheme.  They  applied  a  simple  but 

3  efficient  way  of  avoiding  the  noise  problem  through  neglecting  the  variation  of  the 

4  vertical  diffusion  coefficients  in  the  linearization  of  the  vertical  turbulence  scheme. 

5  Despite  earlier  studies  on  the  parameter  estimation  and  model  verification  (e.g., 

6  Chu  et  al.,  2001),  the  adjoint  model  of  the  turbulence  closure  scheme  has  not  yet  been 

7  thoroughly  investigated  with  either  non-wave  breaking  or  wave  breaking. 

8  Determination  of  wave-affected  parameters  in  the  turbulent  mixing  due  to  breaking 

9  waves  using  the  variation  method  is  selected  as  the  major  objective  of  this  study. 

10  First,  the  upper  layer  temperature  “observations”  are  produced  by  a  “perfect”  model. 

1 1  Second,  a  biased  assimilation  is  conducted  to  identify  the  capability  of  the  variational 

12  method  to  optimally  estimate  the  wave-affected  parameters  in  MY-2.5  turbulence 

13  closure  scheme.  Third,  the  real  temperature  profiles  at  Ocean  Weather  Station  Papa 

14  (OWS  Papa)  are  assimilated  into  the  ocean  model  to  obtain  the  optimal  wave-affected 

15  parameters. 

16  2.  Methodology 

17  2.1  Ocean  boundary  layer  model 

18  Let  (x,  y)  be  the  horizontal  coordinates,  z  the  vertical  coordinate,  and  t  be  the 

19  time.  Following  D'Alessio  et  al  (1998),  equations  governing  the  mean  flow, 

20  temperature,  salinity  in  a  horizontally  homogeneous  ocean  boundary  layer  are  given 
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dv  d  dv 

57\_a/? 

&  H  dz  dz 
&  H  dz 


(1) 


where  w,v  are  the  velocity  components  in  the  x,_y  directions,  respectively,  T  is 
the  potential  temperature,  S  is  the  salinity,  /  is  the  Coriolis  parameter,  Km  and  KH 
are  the  vertical  mixing  coefficients  for  momentum  and  tracers,  respectively. 

The  MY2.5  turbulence  closure  scheme,  widely  used  in  ocean  models  such  as 
POM  and  Regional  Ocean  Modeling  System  (ROMS),  is  a  two-equation  turbulence 
model, 

d<?  ((duV,(dv  2  g  dp 

~wr  ~  +(yr)  >  +  — kh^-~ 

dt  dz  dz  p0  dz 

drl  _  nr,*r  „du. 2  ,  , dv . 2 .  s  „  d 


+ 

| 

rs<  1  Or 

(2) 

q3  W  d  dq2l 

B.l  E,  dz 1  9  dz  ’ 

(3) 

dt  dz  dz  p0  c 

where  q2  is  the  turbulent  kinetic  energy  times  two,  /  is  the  turbulent  macroscale. 

Kq  is  the  vertical  mixing  coefficient  for  turbulence,  p  and  p0  are  the  density  and 

reference  density  respectively, 

W  =  \  +  E2(1/kL)2,  L~x  =  (77  -  z)-'  +  (H  +  z)-' , 
where  k  (=0.41)  is  the  von  Karman  constant,  H  is  the  water  depth,  ij  is  the  free 
surface  elevation,  and  E] ,  E2  and  5,  are  empirical  constants.  The  turbulent  energy 
and  macroscale  equations  are  closed  by 

Km  =  lqSu  ,  K H  =  IqS H ,  Kq  =  IqS q ,  (4) 

where  SM  and  SH  are  the  stability  functions. 

2.2  Wave-affected  parameters 
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1  Wave-affected  parameters  are  included  into  the  surface  boundary  conditions  of 

2  the  two  equation  turbulence  model.  The  first  one  is  the  CB  boundary  condition  for 

3  q 2  (Craig  and  Banner,  1994), 

4  K  =  2aul,  z  =  0  (5) 

9  & 

5  where  uT  is  the  water-side  friction  velocity,  and  a  is  “wave  energy  factor.”  The 

6  second  one  is  for  the  turbulent  macroscale  /  (Terray  et  al.,  1996,  1999), 

7  /  =  max(iczw,/z)  (6) 

8  where  L  is  the  “conventional”  empirical  length  scale,  which  is  calculated 

9  prognostically  by  the  MY2.5  turbulence  closure  scheme;  zwis  the  wave-related 

10  surface  roughness  length,  which  denotes  the  relevant  scale  of  turbulence. 

11  In  the  absence  of  surface  waves,  both  a  and  zw  at  the  surface  are  set  as  zero 

12  in  the  MY2.5  turbulent  closure  scheme  (Blumberg  and  Mellor,  1987).  However,  when 

13  the  effect  of  surface  waves  is  considered,  both  a  and  zw  appear  as  constants  or 

14  vary  with  states  of  surface  waves.  Craig  and  Banner  (1994)  set  a  as  100  for  wave 

15  ages  embracing  very  young  wind  seas  to  fully  developed  situations.  Terray  et  al  (1996) 

16  indicates  that  a  =  150  is  an  adapted  value  under  breaking  waves.  In  thepast,  Kraus 

17  and  Turner  (1967),  Denman  and  Miyake  (1973),  Gaspar  (1988)  also  choose  different 

18  values  of  a  in  their  studies. 

19  Terray  et  al.  (1996),  Burchard  (2001a),  Umlauf  and  Burchard  (2003)  suggest  that 

20  zw  is  the  same  order  as  the  significant  wave  height  (//,).  Further,  Mellor  and 

21  Blumberg  (2004)  summarized  the  work  of  Donelan  (1990),  Smith  et  al.  (1992),  and 
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Janssen  (2001),  and  obtained: 


1  zw  =  /?xl05  X  — ,  (7) 

g 

2  where  g  is  the  gravitational  acceleration,  and  p  is  the  Chamock  parameter  (Chu  and 

3  Cheng  2007),  which  varies  from  P  =  2  (Stacey  1999),  P  =  0.32  (Jones  and  Monismith 

4  2008)  to  p  =  0.56  (Camiel  et  al.  2009)  to  obtain  the  best  performance  in  each 

5  numerical  simulation.  Mellor  and  Blumberg  (2004)  suggested  that  p  ~0(1)  is  deemed 

6  correct  under  breaking  waves.  Smith  et  al.  (1997)  also  indicates  that  P~O(10)  is  too 

7  big  value  to  describe  the  surface  boundary  condition  for  the  turbulent  kinetic  energy. 

8  2.3  Boundary  conditions 

9  The  surface  boundary  conditions  for  q 2  and  /  are  given  by  Eqs.(5)  and  (6). 
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The  bottom  boundary  conditions  of  q 1  and  /  are  given  by 

11 

2  D2/3„  2 

q  —  Bx  uvb 

(8) 

12 

l  -  K20 

(9) 

13 

where  Bl 

=  16.6  (Blumberg  and  Mellor  1987), 

uxb  is 

the  friction  velocity 

14 

associated 

with  the  bottom  frictional  stress.  The 

surface 

and  bottom  boundary 

15 

conditions  of  the  mean  flow  and  tracers  are  represented  by 

16 

17 

18 

cT  Q 
dz  p0Cp 

19 

20 

s=s^ 

dz  dz  A  p0 

>  at 

z  =  0  (10) 

21 

cw  =  (0.75  +  0.067  KDxlO-3  j 
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Cd  =  max( 


x? 
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,0.0025  ) 


>  at 


z  =  -H(x ) 


(11) 


where  Q  is  the  surface  net  heat  flux;  Cp  is  the  specific  heat;  Sobs  is  the 
observation  of  the  sea  surface  salinity;  u,  is  the  friction  velocity  associated  with  the 
wind  stress;  rwx  and  tkv  are  the  x  and  y  components  of  the  wind  stress;  uw  is 
the  wind  speed  at  10  m;  ux  and  u  are  x  and  y  components  of  ul0 ;  xhx  and 
xby  are  the  x  and  y  components  of  the  bottom  frictional  stress;  uh  is  the  bottom 
velocity;  ubx  and  uby  are  the  x  and  y  components  of  uh ;  Cw  and  Cd  are 
drag  coefficients  of  the  wind  stress  and  the  bottom  stress;  and  z0  is  the  bottom 
roughness  parameter,  taken  as  0.01  m. 

2.4  The  variational  analysis 

The  purpose  of  the  variational  analysis  is  to  seek  the  optimal  control  variables  by 
minimizing  a  well-defined  cost  function,  in  which  a  dynamical  model  including  all  the 
control  variables  is  regarded  as  the  strong  constraints  of  the  cost  function.  Within  the 
least-square  framework,  a  general  form  of  the  cost  function  can  be  defined  as 

J(p)  -  W(CX  -  Xobs),{CX  -  Xobs)>  (12) 


where  p  is  the  vector  of  the  control  variables,  X  is  the  solution  of  the  dynamical  model 


dX 

dt 


=  F{X), 
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1  F  is  the  differential  operator.  <  >  is  the  inner  product  in  the  Euclidean  space.  W  is  the 

2  weight  matrix.  Xobs  is  the  observation,  and  C  is  the  projection  operator  from  the 

3  model  space  to  the  observational  space.  Let 

4  J(pobs )  =  min(  p) . 

5  The  optimal  control  variable  pobs  is  obtained  from 

6  VJ(pobs)= 0 

7  with  respect  to  all  control  variables.  Here,  V  is  the  gradient  operator.  The  process  for 

8  the  variational  analysis  can  be  outlined  as  follows: 

9  (a)  Define  a  concrete  cost  function  that  reflects  the  misfit  between  the  control 

10  variables  and  the  available  observations. 

11  (b)  Calculate  the  value  of  the  cost  function  ,J  ( p )  through  integrating  the 
f2  dynamical  model  with  a  fixed  time  step. 

13  (c)  Calculate  the  gradients  of  the  cost  function  with  respect  to  all  control 

14  variables,  VJ(/?). 

1 5  (d)  Minimize  the  cost  function  through  a  minimization  algorithm  according  to 

16  the  value  of  J(p)  and  V/(p). 

17  (e)  Estimate  the  optimal  control  variables  pobs  according  to  the  convergence 

1 8  criterion  of  the  process  of  the  minimization. 

19  For  executing  the  above  process  of  the  variational  analysis  (a)-(e),  V/(/?)  should  be 

20  obtained  in  advance,  and  in  general, it  is  calculated  by  the  adjoint  model  of  the 

21  linearized  dynamical  model.  To  the  first  order  the  Taylor  expansion  of  J(p)  is 

22  given  by 

23  J(p)  =  J(p0)  +  SJ(p)  (13) 

24  where  3J(p)  is  the  variation  of  the  J(p)  .On  the  one  hand,  3J(p)  is  given  by  the 


li 


1  definition  of  the  variation: 

2  S/(p)=l<VM,X>  (14) 

3  On  the  other  hand,  SJ{p )  can  also  be  written  according  to  the  Eq.  (12) 

4  8J(p)  =  i  £<  W^SX,  CX  -  Xobs  >  +  i  [<  W  {CX  -  Xobs  ),^-8X>.(  15) 

5  With  the  symmetry  of  the  inner  product  as  well  as  a  constant  W  matrix,  Eq.  (15)  can 

6  be  rewritten  as 

7  aj{p)  =  \T<W{CX-Xobs),^8X>  (16) 

J0  oX 

8  Let  - =A(X),  thus  Eq.  (15)  can  be  given  as 

dX 

9  &J(p)  =  j<  W(CX  -  Xobs),  A(X)8X  >  (17) 

10  where  A(X)  is  called  as  the  tangent  linear  operator.  Eq.(17)  can  be  transposed 

1 1  according  to  the  definition  of  the  adjoint  operator: 

12  8J(p)  =  £  <  WA*(X)(CX  -  Xobs),SX  >  (18) 

13  where  A*(X)  is  called  as  the  adjoint  operator  of  A(X) .  Compared  with  Eq.  (14), 

14  VxJ(p0)  can  be  described  by: 

15  VxJ(p0)  =  WA\X)(CX-Xobs).  (19) 

16  According  to  Eq.  (19),  the  gradient  of  the  cost  function  with  respect  to  the  control 

17  variables  can  be  calculated  using  the  adjoint  model.  The  difference  CX  -Xobs  is 

18  regarded  as  an  external  forcing  of  the  adjoint  model. 

19  2.5  The  adjoint  model 

20  The  dynamical  model  composed  by  Eq.(  1) — (1 1)  can  be  summarized  in  a  general 
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form  as 


2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 


dx 
dt 

x\  =  xn 

K  0 

x(t) \r  =  y(t) 


(20) 


where  x  is  the  vector  of  model  state  variables,  including  u,v,T,  S,  q 2  and  q2l ; 
x0  is  the  model  states  at  initial  time  t0 ;  ,  and  y{t)  is  the  boundary  condition  on  r 

respectively. 

The  tangent  linear  model  of  Eq.  (20)  can  be  written  by 

dx'  _  dF(x)  , 
dt  dx 

X'L=X0  (21) 


x\t)\r  =  y'(t) 

where  the  prime  is  the  perturbations  of  the  state  variables. 

For  the  two  vectors  w  and  z  in  the  Euclidean  space,  the  adjoint  operator  L*  of  the 
linear  operator  L  can  be  defined  as: 


* 

<z,  Lw>=<L  z,  w> 

In  the  Euclidean  space,  L  is  the  transpose  of  L,  namely  L*=  LJ.  The  adjoint  model 
corresponding  to  (21)  is  given  by 


dx  _  dF(x)  T 
dt  dx 


x  =0 


x 


«t  =  o 


(22) 


where  x  represents  the  adjoint  variables  and  tE  is  the  end  time  in  the  temporal 


integration  of  Eq.  (20).  The  negative  sign  in  the  right  side  of  the  first  equation  in  (22) 
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1  indicates  that  the  adjoint  model  integrates  backward  in  time.  When  the  adjoint  model 

2  integrates  backward  to  the  initial  time  t0,  the  corresponding  3c  | t_t  is  the  gradient  of 

3  the  cost  function  with  respect  to  the  state  variables. 

4  The  discretized  adjoint  model  that  computes  the  gradient  of  the  cost  function  can 

5  be  developed  directly  from  the  discretized  dynamical  model  including  Eqs.(l)~(ll). 

6  In  practical  application,  the  source  code  of  the  adjoint  model  is  constructed  by 

7  combining  the  Tangent  and  Adjoint  Model  Compiler  (TAMC)  developed  by  Giering 

8  and  Kaminski  (1998)  and  a  hand-coding  correction.  First,  the  adjoint  code  is 

9  generated  by  TAMC  to  avoid  man-made  errors  and  negligence,  which  are  extremely 

10  easy  to  happen  during  the  direct  coding.  Second,  hand-coding  correction  is  conducted 

11  to  correct  the  AMC-generated  code  and  control  the  adjoint  code  structure.  The  errors 

12  in  the  adjoint  code,  which  are  induced  from  some  irregular  expressions  of  the  forward 

13  numerical  model  such  as  the  partial  array  assignment  and  iterative  use  of  intermediate 

14  arrays,  are  corrected  through  the  hand  coding.  Finally,  through  the  hand-coding 

15  correction,  values  of  many  intermediate  results  in  the  adjoint  model  are  recorded  into 

16  memory  instead  of  recomputation  to  shorten  run  time  of  the  adjoint  model,  and  some 

17  local  variables  and  arrays  are  transferred  to  global  attribute  to  improve  the  run 

18  efficiency  of  the  adjoint  model. 

19  Once  the  cost  function  and  its  gradient  are  obtained  from  the  dynamical  model 

20  and  associated  adjoint  model,  the  minimization  process  is  implemented  to  minimize 

21  the  cost  function  through  iterating  the  values  of  the  control  variables  (T"  ,T"~\  a 

22  and  J3  )  with  the  limited  memory  Broyden-Fletcher-Glodfarb-Shanno  (BFGS) 
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quasi-Newton  minimization  algorithm  (Liu  and  Nocedal,  1989).  During  the 
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minimization  process,  the  maximum  of  a  is  set  to  1000,  and  the  maximum  of  /?  is 
set  to  10  according  to  Mellor  and  Blumberg  (2004)  and  Smith  et  al.  (1997).  The 
minima  of  the  two  wave-affected  parameters  are  set  to  zero  to  keep  realistic 
physical  conditions.  The  minimization  process  is  repeated  until  the  convergence 
criterion  of  the  gradient  is  reached.  At  that  time,  the  optimal  values  of  the  control 

variables  are  obtained. 

2.6  Cost  function 

The  cost  function  is  defined  by 


J(Tn,Tn~\a,P)  =  -{T" -TbnY  Bp  (Tn  -T£)  +  —  (T  -7),  )  Bp(T  -  7)  ) 

1  M  N 

+  ^TJTJVjM’P)-ToJTR-\T].{a,P)-Tobs) 

^  j= 1  *=1 


(23) 


where  the  first  two  terms  in  the  right  side  represent  the  background  error  terms  that 
measure  the  misfit  between  the  model's  initial  field  and  the  background  field.  T" 


and  T"  1  are  the  initial  temperature  values  at  the  nth  and  (n-l)m  time  step  respectively, 


\th  , 


which  will  be  estimated  optimally  via  the  variational  method.  7)”  and  7)"  1  are  the 


background  temperature  values  at  the  nth  and  (n-l)th  time  steps  respectively,  which  can 
be  derived  from  the  model  run.  Both  temperatures  at  the  two  consecutive  time  steps 
are  considered  as  the  control  variables  due  to  the  utilization  of  the  leapfrog  time 
differencing  scheme  with  the  Asselin-Robert  time  filter  (Robert,  1966).  Otherwise, 
initial  shocks  of  the  model  states  are  likely  to  be  produced  during  the  variational 
estimation  because  of  the  inconsistence  of  the  initial  values  at  the  two  time  steps.  Bi 


and  Bt  are  the  error  covariance  for  T"  and  T"  1  respectively,  for  simplicity,  both  B i 
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and  i?2use  diagonal  matrices,  whose  values  of  the  diagonal  components  are  set  to  10’4 
in  this  study.  The  third  term  denotes  the  observation  of  the  temperature  at  certain  time 
intervals  within  the  assimilation  window,  where  T, ,  and  Tobs  are  the  simulated  and 
observed  temperature  at  location  i  and  time  level  j.  N  and  M  are  the  number  of  grid 
points  over  the  ocean  and  the  number  of  time  levels  of  observations.  R  is  the  error 
covariance  for  the  observations,  which  also  uses  the  same  diagonal  matrix  as  that  of 
Bi. 

Wave-affected  parameters  a  and  /?  are  expressed  implicitly  in  Eq.  (23),  which  are 
regarded  as  the  independent  variables  of  77 , .  Therefore,  the  value  of  the  cost  function 
can  be  obtained  when  the  model  integrates  for  n  time  steps  with  the  known  initial 
values  of  T"  ,Tn~l ,  a  and  f3 .  The  cost  function  has  the  following  form  if  the 
wave-affected  parameters  a  and  /?  have  background  values  (a*,  /?/,), 


J(Tn  ,Tn\a,P)  =  ^{Tn -T”f  B;l(Tn  -  Tb" )  +  ^  (T"'1  -  T^'  f  B2 1  (r  1  -  T/"'1 ) 


1  M .  J1  1  1 

+  -YJYJ(TjM,f3)-T0JTR-\TjAa,f])-Tobs)+2KM-ab)2  +  -Kfi(j3-j3b) 


(24) 


where  Ka  and  K f,  are  coefficients  controlling  the  best  fits  for  data.  In  this  study, 
we  use  the  first  form  of  the  cost  function  (23)  for  avoiding  the  complexity  of  the  cost 
function. 

3.  Synthetic  experiments 

3.1.  “Truth”  model  simulation 

Table  1  lists  all  the  assimilation  experiments  and  model  simulations  within  an 
identical  synthetic  experiment  framework.  The  “truth”  model  consists  of  Eq.(l)-(3) 
with  a  =  200  and  /?=2.  All  the  6  equations  from  (1)  to  (3)  are  discretized  using  the 
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same  implicit  method  as  POM.  The  maximum  depth  is  set  to  250  m,  with  60  vertical 
levels.  The  first  20  vertical  levels  are  0.0,  0.5,  1.0,  1.5,  2.0,  4.0,  6.0,  8.0,  10.0,  12.0, 
14.0,  16.0,  18.0,  20.0,  22.0,  24.0,  26.0,  28.0,  30.0,  35.0  m.  The  time  step  is  1-hour. 
The  model  initial  state  is  from  Jan.  1,  1961,  including  temperature  and  salinity, 
derived  from  the  real  observation  at  OWS  Papa.  The  model  is  forced  by  the 
observational  10-minutes  momentum  and  heat  fluxes  from 
http://www.mnel.noaa.gov/stnp/data.html. 

Starting  from  the  initial  conditions  (Jan.  1,  1961),  the  "truth"  model  is  run  for  6-yr 
to  generate  time  series  of  the  "truth"  with  the  first  5-yr  as  the  spin-up  period.  The  time 
of  the  "observations"  of  T  is  from  Aug.  1,  1966  to  Aug.  30,  1966.  The 
"observations"  of  T  are  produced  through  sampling  the  "truth"  states  at  1-hour 
observational  frequencies.  The  "observation"  locations  of  T  are  consistent  with  those 
of  the  model  vertical  grids. 

3.2.  Biased  simulation 

The  biased  simulation  uses  the  same  “truth”  model,  but  with  different  parameter 
settings.  Therefore,  the  difference  between  the  biased  simulation  and  “truth”  model 
leads  to  the  effect  of  the  “incorrect”  parameter  settings.  Fig.l  shows  the  simulated 
daily  temperature  at  OWS  Papa  in  1966.  The  sea  surface  temperature  (SST)  from  the 
biased  simulation  with  (a,j3)  =  (100,  1)  is  higher  than  that  by  the  “truth”  model 
simulation  with  (a  ,J3)  =  (200,  2),  and  the  maximum  difference  of  the  SST  between 
the  two  simulations  occurs  in  summer,  namely  from  the  200th  day  to  the  240th  day 
(solid  line  vs.  dash  line  in  Fig. la).  Obvious  difference  of  the  temperature  at  10  m 
depth  in  the  two  simulations  also  remains  (Fig. lb).  The  wave-affected  parameters  are 
half  smaller  in  the  biased  simulation  than  in  the  “truth”  model  simulation,  which 
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1  suggests  that  the  turbulent  kinetic  energy  is  too  weak  to  mix  the  surface  and 

2  subsurface  water  well  in  the  biased  simulation.  After  the  240th  day  (fall  and  winter), 

3  the  temperature  decreases  gradually  due  to  the  convective  mixing  induced  by  the 

4  surface  cooling.  The  temperatures  at  the  surface  and  10  m  depth  in  the  biased 

5  simulation  remain  higher  than  the  counterpart  in  the  “truth”  model  simulation  due  to 

6  the  insufficient  wave-enhanced  mixing  in  the  biased  simulation.  Below  20  m,  the 

7  effect  of  the  wave-affected  parameters  on  the  temperature  is  not  evident  in  summer 

8  (solid  line  vs.  dash  line  in  Figs.lc  and  Id),  which  indicates  that  the  turbulent  kinetic 

9  energy  generated  by  the  breaking  surface  gravity  waves  is  dissipated  only  near  the  sea 

10  surface  and  does  not  penetrate  into  the  deeper  waters.  The  maximum  difference  in 

11  temperature  at  30  m  from  the  two  simulations  occurs  in  the  fall  (after  the  250th  day) 

12  with  temperature  higher  in  the  biased  simulation  than  in  the  “truth”  model  simulation. 

13  Although  the  wave-affected  parameters  do  not  directly  affect  the  temperature  in  the 

14  deeper  layers  in  summer,  it  can  affect  the  temperature  indirectly  by  the  SST  due  to  the 

15  subsequent  convective  cooling  in  autumn  and  winter.  Thus,  the  wave-affected 

16  parameters  directly  impact  on  the  temperature  near  the  sea  surface  in  summer,  and 

1 7  indirectly  impact  on  the  temperature  in  the  deeper  layers  in  autumn  and  winter. 

18  We  intend  to  investigate  if  the  wave-affected  parameters  in  a  two-equation 

19  turbulence  model  can  be  estimated  effectively  through  assimilating  the  temperature 

20  data  into  an  ocean  boundary  layer  model  with  the  variational  method.  In  addition,  we 

21  want  to  understand  how  well  the  model  state  estimation/forecast  can  be  improved 

22  through  the  estimated  wave-affected  parameters.  In  the  next  subsection,  a  series  of 
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1  synthetic  experiments  are  carried  out  to  address  the  issues. 

2  3.3.  Parameter  estimation 

3  Fig.2  shows  a  flowchart  of  the  wave-affected  parameter  estimation  with  the 

4  variational  method.  The  process  for  the  wave-affected  parameter  estimation  is 

5  outlined  as  follows: 

6  (a)  Begin  with  the  initial  field  on  Aug.  1,  1966  and  use  the  different  values  of 

7  wave-affected  parameters  from  the  “truth”  for  the  biased  simulation. 

8  (b)  Integrate  the  model  Eqs.  (1)  -(3)  forward  to  a  fixed  time  window  ATw  and 

9  calculate  the  value  of  the  cost  function  J(T',T'  \cx,j3)  using  Eq.  (23). 

10  (c)  Integrate  the  adjoint  model  backward  in  time  and  calculate  the  values  of  the 

1 1  gradient  of  the  cost  function  with  respect  to  the  control  variables  VJ. 

12  (d)  With  the  values  of  the  cost  function  J{T ,T~\oL,fi)  and  the  gradient  VJ,  use 

13  the  BFGS  algorithm  to  obtain  the  new  values  of  the  control  variables, 

14  namely,  the  two  wave-affected  parameters  a ,  (5  and  initial  upper  layer 

15  temperature  fields  T"  ,T"~l . 

16  (e)  With  the  updated  control  variables  from  process  (d),  repeat  processes  (b)  -  (d) 

17  until  the  convergence  criterion  for  the  minimization  is  satisfied.  The 

1 8  convergence  criterion  is  defined  as 

19  |W|/|W0|<0.01. 

20  The  solution  of  the  control  variables  that  satisfies  the  convergence  criterion  is 

2 1  regarded  as  the  optimal 

22  solution. 


19 


1  (f)  Integrate  the  model  Eqs.  (1)  -(3)  to  the  fixed  time  window  ATw  using  the 

2  optimal  solution  derived  from  process  (e),  and  results  are  regarded  as  the  new 

3  initial  fields  for  the  next  integration. 

4  (g)  Use  the  new  initial  fields  derived  from  the  process  (f)  and  the  optimal 

5  wave-affected  parameters  derived  from  process  (e),  iterate  the  processes  (b) 

6  to  (f)  to  obtain  time  series  of  wave-affected  parameters  a  and  f5 . 

7  Fig.  3  shows  the  time  series  of  a  and  /?  during  the  parameter  estimation  (PE) 

8  described  in  Table  1,  where  both  assimilation  window  and  frequency  are  set  to  24 

9  hours  and  the  assimilation  depth  is  set  to  30  m.  Therefore,  the  processes  (b)~(g)  are 

10  executed  30  times  to  obtain  time  series  of  a  and  /?  from  Aug.  1,  1966  to  Aug. 

1 1  30,  1966.  Fig.3b  shows  that  /?  converges  to  its  “truth”  value  (dash  line)  after  9  days, 

12  while  a  converges  to  its  “truth”  value  (dash  line  in  Fig.3a)  after  about  15  days. 

13  Results  show  the  wave-affected  parameters  in  the  high  order  turbulent  model  can  be 

14  estimated  successfully  using  the  upper  layer  temperature  observations  through  the 

15  variational  control  technique.  For  each  cycle  of  the  parameter  estimation  in  the  30 

16  days,  the  process  of  the  minimization  is  iterated  until  the  convergence  criterion  of  the 

17  gradient  is  satisfied.  Fig. 4  shows  the  dependence  of  the  cost  function  and  the  norm  of 

18  the  gradient  on  the  number  of  iterations  on  Aug.  2,  1966.  The  value  of  the  cost 

19  function  decreases  rapidly  from  4.3  to  0.8  within  first  5  iterations,  and  keeps  the  low 

20  value  (0.8)  steadily  after  the  5th  iteration  (Fig.  4a).  However,  the  norm  of  the  gradient 

21  oscillates  dramatically  to  search  the  optimal  declining  direction  of  the  gradients.  The 

22  norm  of  the  gradient  goes  stable  after  the  130th  iteration  (Fig.  4b).  The  minimization 
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1  process  stops  after  180th  iterations,  indicating  the  local  minima  of  the  wave-affected 

2  parameters  for  that  day. 

3  Fig.  5  is  the  temporal  variations  of  the  natural  logarithm  of  the  cost  function  at 

4  OWS  Papa  from  Aug.  1  to  Aug.  30,  1966.  The  cost  function  (red  line)  decreases 

5  dramatically  in  the  first  5  days,  then  decreases  gently  in  the  following  25  days.  Both 

6  the  background  term  (blue  line  in  the  Fig.5)  and  the  observation  term  (black  line  in 

7  the  Fig.5)  of  the  cost  function  have  a  similar  pattern  with  the  total  cost  function.  The 

8  two  terms  almost  converge  to  the  same  value  after  the  10th  days,  indicating  the 

9  estimated  initial  temperature  fields  reach  a  balance  between  the  background 

10  temperature  and  the  observation. 

1 1  The  temporally  varying  wave-effected  parameters  (a  ,/3)  estimated  from  their 

12  different  initial  values  on  Aug.  1,  1966  (Fig.  6)converge  to  their  “truth”  values  within 

13  one  month  through  the  parameter  optimization  with  the  variational  approach.  It 

14  clearly  shows  that  the  variational  assimilation  approach  is  feasible  for  the 

15  wave-affected  parameter  optimization  with  different  initial  parameter  values. 

16  To  evaluate  the  effect  of  the  noise  in  the  temperature  observation  on  the 

17  wave-affected  parameter  estimation,  based  on  the  PE  experiment,  the  white  noises 

18  with  different  standard  deviation  are  added  to  the  temperature  observation.  Table  2 

19  shows  the  dependence  of  the  optimally-estimated  (a,j3)  on  the  error  standard 

20  deviation  of  temperature  observation.  The  relative  error  of  optimally-estimated  a 

21  decreases  from  96.9%  to  60%,  and  the  relative  error  of  optimally-estimated  J3 

22  decreases  from  99.1%  to  94.3%  as  the  error  standard  deviation  in  temperature 
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1  observation  increases  from  0.001  to  0.05K.  It  implies  that  the  effect  of  observational 

2  noise  on  the  estimation  is  more  severe  on  a  is  than  on  [i ,  which  means  that  it  is 

3  more  difficult  to  pick  up  the  useful  signal  when  the  noise  dominates  the  cost  function 

4  and  corresponding  gradients  during  the  parameter  estimation  of  a  .  When  the 

5  standard  deviation  of  temperature  observation  increases  to  0.5K,  both  relative  errors 

6  of  the  optimally-estimated  a  and  j5  are  below  50%,  which  indicates  that  the  level 

7  of  the  noise  is  not  acceptable  for  assimilation  purposes. 

8  To  explore  if  the  wave-affected  parameters  can  be  estimated  correctly  only  using 

9  the  SST  data,  the  second  assimilation  experiment,  PE_SST,  is  conducted,  in  which 

10  only  the  SST  observations  are  assimilated  into  the  biased  simulation  model.  Neither 

11  a  (Fig.  7a)  nor  /?  (Fig.  7b)  reaches  their  “truth”  values  (dashed  curve)  due  to  the 

12  poor  constraint  of  the  observation.  When  only  the  SST  observations  are  assimilated, 

13  the  subsurface  temperature  cannot  be  estimated  accurately.  Under  this  condition,  the 

1 4  two  parameters  will  be  adjusted  to  the  optimal  values  to  fit  the  inaccurate  temperature 

15  values  to  the  greatest  extent  within  a  fixed  time  window,  rather  than  to  converge  to 

16  “truth”  values.  Therefore,  the  subsurface  temperature  observations  are  essential  for 

17  estimating  a  and  /?  reasonably  well. 

18  The  dependence  of  the  optimally-estimated  a  (Fig.  8a)  and  /?( Fig.  8b)  on  the 

19  assimilation  window  and  frequency  is  investigated  using  different  values  from 

20  Aug.l  to  Aug.  30,  1966  (Fig.  8).  When  the  assimilation  window  and  frequency  are  48 

21  hours  and  72  hours,  both  parameters  converge  to  their  respective  “truth”  values  (see 

22  black  and  blue  lines  in  Figs.  8a  and  8b).  Flowever,  when  the  assimilation  window  and 
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1  frequency  reach  96  hours  and  120  hours,  neither  a  nor  [3  converges  to  their 

2  “truth”  values  within  one  month,  which  can  be  seen  from  the  red  and  pink  lines  in 

3  Figs.  8a  and  8b.  It  clearly  shows  that  the  parameter  updating  with  the  observation  can 

4  improve  the  state  estimation  of  the  next  cycle,  and  the  improved  state  estimation 

5  further  enhances  the  quality  of  parameter  estimation  for  the  next  cycle  of  parameter 

6  correction.  When  the  assimilation  window  and  frequency  are  set  to  120  hours,  the 

7  state-parameter  optimization  is  performed  only  in  6  cycles  within  one  month. 

8  Although  the  cost  function  decreases  gradually,  which  can  be  seen  from  the  dash 

9  curve  in  Fig.  9,  the  control  variables  (the  initial  temperature  7a nd  two  parameters  a,  (:i) 

10  are  not  estimated  reasonably  well.  In  contrast,  when  the  assimilation  window  and 

11  frequency  are  set  to  24  hours,  just  as  the  PE  experiment,  the  state-parameter 

12  optimization  can  be  performed  in  30  cycles  within  one  month,  and  the  cost  function 

13  can  reach  quasi-equilibrium  after  10  days  (solid  curve  in  Fig. 9). 

14  The  incorrect  convergence  of  (a,  P)  suggests  that  the  initial  temperature  field  is 

15  not  adjusted  well  enough,  which  is  regarded  as  the  source  of  noise  during  parameter 

16  estimation  using  the  variational  method.  Therefore,  it  is  hard  to  obtain  the  accurate 

17  values  of  (a,  p)  before  the  state  variables  (  Tn  and  T"~l)  attain  the  adequate  accuracy. 

18  For  better  understanding  the  issue,  two  other  experiments  are  carried  out,  in  which  /3 

19  is  regarded  as  the  only  control  variable.  The  experiment  PE_/?_TI  described  in  Table 

20  1  uses  the  “perfect”  initial  field  that  is  generated  by  the  “truth”  model  with  the  “truth” 

21  values  of  a  and  /?,  the  other  experiment,  PE_/?_BI,  uses  the  “biased”  initial  field 

22  that  is  generated  by  the  biased  simulation  with  the  “biased”  values  of  a  and  (3 . 
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1  Table  3  shows  the  evolution  of  the  cost  function,  norm  of  the  projected  gradient  and 

2  value  of  /?  with  respect  to  the  number  of  iterations  in  PE/iTI.  The  parameter  [i 

3  reaches  its  “truth”  value  at  the  3rd  iteration.  The  convergence  criterion  of  the  gradient 

4  is  satisfied  at  the  4th  iteration.  However,  (3  estimated  from  PE_/?_BI  cannot 

5  converge  to  its  “truth”  value  (Table  4).  After  the  convergence  criterion  of  the  gradient 

6  is  satisfied  at  the  6th  iteration,  f3  reaches  3.302335,  which  is  different  from  the 

7  “truth”  value  2.0.  Although  f3  from  PE_/?_BI  cannot  converge  to  its  “truth”  value,  it 

8  reaches  its  optimal  value  to  compensate  the  error  derived  from  the  “biased”  initial 

9  filed  during  minimizing  the  model-observation  misfit. 

10  In  fact,  in  a  3D  ocean  circulation  model,  model  biases  arise  from  the  imperfect 

1 1  dynamical  core  and  empirical  physical  schemes  even  if  the  initial  field  is  perfect.  With 

12  a  biased  initial  field  alone,  one  expects  that  the  parameter  optimization  can 

13  compensate  both  numerical  and  physical  deficiencies  of  numerical  model  and  enhance 

14  the  performance  of  the  model  simulation  to  certain  degree.  Under  this  situation, 

15  parameters  can  only  converge  to  their  optimal  value,  instead  of  the  “truth”  values.  In 

16  the  next  section,  real  temperature  profiles  from  OWS  Station  Papa  will  be  assimilated 

17  into  the  assimilation  model  to  obtain  the  optimal  wave-affected  parameters  (a,  p). 

18  4.  Real  experiment 

19  The  Papa  Station  locates  in  the  North  Pacific  at  (145°W,  50°N),  where  the 

20  currents  are  relatively  weak  and  the  local  mixing  modulates  mainly  the  dynamical 

2 1  process  in  the  upper  ocean  in  summer.  The  observed  temperature  profiles  from  Aug 

22  1  to  Aug  31,  1966  at  the  site  have  3h  time  interval  and  a  coarser  vertical  resolution 
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1  (5  m)  than  the  model  grid  points.  There  are  7  observational  layers  totally  in  the  upper 

2  30  m,  namely  0,5,10,15,20,  25,30  m.  Linear  interpolation  is  used  to  fill  the  spatial  gap 

3  between  the  modeled  and  the  observational  data. 

4  Table  5  lists  all  the  assimilation  experiments  and  model  simulations  within  the 

5  real  experiment  framework.  First,  a  control  run  without  assimilating  any  observational 

6  data,  is  called  control  (CTRL)  to  serve  as  the  reference  for  the  evaluation  of 

7  assimilation  experiments.  The  initial  temperature  and  salinity  are  taken  from  those  at 

8  00:00  GMT  on  Jan  1,  1961,  and  linearly  interpolated  to  model  grids.  The  high 

9  resolution  (1/6°)  surface  observed  data  (momentum  and  net  heat  fluxes)  at  the  site  are 

10  used  to  force  the  model.  Fig.  10a  shows  the  daily  observed  (red  curve)  and  simulated 

1 1  sea  surface  temperature  from  CTRL  (black  dashed  curve)  at  the  OWS  Papa  on  Aug., 

12  1966.  The  simulated  SST  is  higher  than  the  observed  SST  by  about  3°C  (black  dashed 

13  curve  vs.  red  curve).  At  the  same  time,  the  simulated  mixed  layer  depth  from  CTRL  is 

14  shallower  than  the  observation  by  more  than  10  m  (black  dashed  curve  vs.  red  curve 

15  in  Fig.  10b).  The  optimal  values  of  (a,  p)  are  estimated  with  the  variational  method  to 

16  mitigate  the  bias  between  the  model  and  the  observation  using  the  real  summer 

1 7  temperature  data. 

18  The  real  parameter  estimation  (RPE)  is  described  in  the  second  row  of  Table  5. 

19  The  initial  field  is  generated  from  the  results  on  Aug.  1,  1966  simulated  by  the 

20  “truth”  model  in  the  above  synthetic  experiments.  The  initial  values  of  (a,  P)  are  also 

21  consistent  with  those  in  the  “truth”  model  simulation.  The  length  of  the  assimilation 

22  window  is  set  to  3  days  (8  real  observational  temperature  profiles  in  each  day,  totally 

23  24  profiles  within  3  days)  and  the  assimilation  depth  is  30  m.  The  process  of  PE  is 
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1  similar  to  the  process  described  in  Section  3,  but  with  the  real  temperature 

2  observations  at  OWS  Papa  in  Aug.,  1966.  Table  6  shows  the  evolution  of  the  cost 

3  function,  a  and  f3  with  respect  to  the  number  of  iteration  for  RPE.  After  the  8th 

4  iteration,  the  normalized  cost  function  decreases  to  5%  of  its  initial  value.  The  optimal 

5  values  of  a  and  f3  reach  107.48  and  3.98  respectively.  The  SST  from  RPE  has  a 

6  significant  improvement  compared  to  the  simulated  SST  from  CTRL  (black  solid 

7  curve  vs.  black  dashed  curve  in  Fig.  10a),  whose  values  are  basically  consistent  with 

8  those  of  the  observations  (black  solid  curve  vs.  red  curve  in  Fig.  10a).  The  mixed 

9  layer  depth  is  also  more  accurate  from  RPE  than  from  CTRL  (Fig.  10b).  However, 

10  some  discrepancy  in  the  mixed  layer  depth  still  exist  between  RPE  and  the 

11  observation.  This  is  because  too  many  factors  modulate  the  complicated 

12  thermodynamic  processes  of  the  upper  mixed  layer  besides  the  surface  gravity  waves, 

13  such  as  horizontal  advection,  internal  waves,  upwelling,  and  entrainment.  Many 

14  physical  processes  are  not  enclosed  in  the  simple  ocean  boundary  layer  model.  The 

15  optimal  values  of  the  parameters  can  only  compensate  some  model  bias,  but  not  all. 

16  However,  the  result  from  RPE  indicates  that  the  variational  estimation  of 

17  wave-affected  parameters  can  indeed  reduce  model  biases  and  improve  the  model 

1 8  capability  in  the  upper  ocean. 

19  To  explore  the  impact  of  parameter  estimation  on  model  simulation,  two 

20  validation  experiments,  RSEPo  and  RSEPd,  are  conducted.  The  “optimal” 

21  parameters  estimated  from  RPE  are  used  in  RSE  Po,  and  the  default  values  of  the 

22  parameters  from  CTRL  are  used  in  RSE  Pd.  In  addition,  both  experiments  use  the 

23  same  initial  fields  on  Aug  31,  1966,  which  are  derived  from  RPE.  Fig.  11  shows  the 

24  observed  (red  curve)  and  simulated  SST  from  RSE  Po  (black  solid  curve)  and 
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1  RSEPd  (black  dash  curve)  at  OWS  Papa  from  Sept.  1  to  Sept  30,  1966.  The 

2  simulated  SST  is  more  consistent  with  the  observations  from  RSEPo  than  from 

3  RSE  Pd.  The  simulated  twice  monthly-averaged  turbulent  kinetic  energy  q 2  (Fig. 

4  12a),  and  vertical  mixing  coefficient  for  temperature  KH  (Fig.  12b)  at  OWS  Papa  in 

5  Sept  1966  are  much  larger  for  all  depths  in  RSE  Po  (solid  curve)  than  in  RSE  Pd 

6  (dashed  curve).  The  enhanced  KH  in  the  upper  30m  depth  in  RSE_Po,  due  to  the 

7  improvement  of  the  turbulent  kinetic  energy  calculation  ,  mixes  the  momentum  from 

8  the  winds  downward  through  the  water  column  and  makes  it  more  vertically 

9  homogeneous.  It  indicates  that  the  model  performance  can  be  effectively  improved 

10  using  the  optimal  parameters.  However,  more  accurate  model  simulations  are  needed 

1 1  using  the  optimal  values  of  parameters  via  the  variational  methods  repeatedly  at  the 

12  certain  time  intervals  with  more  available  observations. 

13  5.  Discussion  and  conclusion 

14  The  upper  layer  temperature  data  is  assimilated  into  an  ocean  surface  boundary 

15  layer  model  to  estimate  the  wave-affected  parameters  (a,  p)  employed  in  the  MY2.5 

16  two-equation  turbulence  model  using  the  variational  method.  Within  an  identical 

17  synthetic  experiment  framework,  the  “truth”  values  of  the  wave-affected  parameters 

18  in  the  high  order  turbulence  model  can  be  retrieved  successfully  when  the  assimilation 

19  window,  the  assimilation  frequency,  and  the  assimilation  depth  are  set  appropriately. 

20  The  observational  temperature  profiles  at  the  OWS  Station  Papa  are  also  assimilated 

21  to  correct  the  model  bias  arisen  from  multiple  sources.  By  fitting  the  model  results  to 

22  the  observations  using  the  variational  method,  the  optimal  temperature  field  can  be 
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1  obtained  in  the  upper  30  m  through  adjusting  the  wave-affected  parameters  to  their 

2  optimal  values.  Wave-affected  parameters  estimation  using  the  variational  method  can 

3  compensate  in  part  the  numerical  and  physical  deficiencies  of  the  model  in  the  upper 

4  ocean.  However,  It  should  also  be  noted  that  the  optimal  values  of  the  wave-affected 

5  parameters  are  not  the  so  called  “truth”  values.  The  “optimal”  values  of  the 

6  wave-affected  parameters  in  real  applications  are  only  applicable  to  the  specific 

7  time  period,  location,  and  model.  The  “optimal”  values  should  vary  temporally  and 

8  spatially  rather  than  being  constants,  which  can  be  obtained  by  using  the  variational 

9  methods  repeatedly  at  the  certain  time  intervals  and  available  observations  (Peng  et  al, 

10  2012).  Although  the  optimal  values  of  the  wave-affected  parameters  are 

11  model-dependent  (initial  fields,  time  window  of  assimilation,  model  configuration, 

12  etc.),  they  can  indeed  mitigate  the  model  biases  from  multiple  sources,  and  obviously 

13  improve  the  performance  of  the  model  simulation.  Besides  the  wave  breaking 

1 4  parameters,  other  parameters  in  the  wave-related  processes  can  also  be  introduced  into 

15  the  model  (which  is  compatible  with  those  pertinent  to  the  wave  breaking  )  to 

16  estimate  their  optimal  values. 

17  In  general,  the  complex  turbulent  closure  models  are  empirical  and  full  of 

18  uncertainty  in  an  ocean  circulation  model.  Due  to  the  high  nonlinearity  and 

19  discontinuity  of  the  vertical  turbulence,  it  is  more  difficult  to  linearize  the  complicated 

20  turbulence  closure  scheme  than  to  linearize  the  momentum  and  tracer  equations. 

21  Wave-affected  parameters  in  high  order  turbulence  closure  schemes  can  modulate 

22  distinctly  the  vertical  structure  in  the  upper  ocean.  Therefore,  it  is  essential  to  estimate 
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1  their  optimal  values  using  observations  deployed  in  the  upper  ocean  through  some 

2  robust  data  assimilation  methods  such  as  the  variational  method  or  the  ensemble 

3  Kalman  filter.  Now,  satellite  remote  sensed  SST  data  and  in-situ  temperature  data 

4  (such  as  the  Argo  floats)  can  provide  a  mass  of  temperature  observations  in  upper 

5  oceans.  Therefore,  the  optimal  geographic-dependent  distribution  of  the  wave-affected 

6  parameters  in  a  high  order  turbulence  closure  scheme  can  be  obtained  using  the  4DVar 

7  that  assimilates  the  upper  layer  available  temperature  data  into  ocean  circulation 

8  models. 

9  Appendix  A.  Sensitivity  of  simulated  temperature  to  parameters 

10  It  is  essential  to  investigate  model  sensitivities  with  respect  to  parameters  being 

11  estimated  before  parameter  estimation.  Fig.  A1  shows  the  dependence  of  the  cost 

12  function  on  a  and  j5 .  It  increases  with  the  increasing  a  and  [i  in  general. 

13  However,  the  local  minimum  of  the  cost  function  can  be  found  near  the  region  in 

14  which  both  a  and  (5  reach  their  default  values  (see  Fig.  Alb).  The  existence  of  the 

15  local  minimum  indicates  that  it  is  likely  to  estimate  the  optimal  values  of  a  and  J3 

16  if  the  values  of  the  gradient  with  respect  to  the  parameters  can  be  calculated  correctly 

17  in  all  the  numerical  iterations  by  the  adjoint  model. 

18  The  ensemble  spread  of  T  is  used  to  evaluate  the  relevant  sensitivities 

19  quantitatively.  For  a  and  J3 ,  100  Gaussian  random  numbers  are  generated  with  the 

20  standard  deviation  being  5%  of  the  default  value  and  superimposed  into  the  parameter 

21  being  perturbed,  while  the  other  parameter  remain  unperturbed.  All  the  100  ensemble 

22  members  are  started  from  the  same  initial  conditions  (Jan.  1,  1961).  The 

23  biased-simulation  model  is  integrated  up  to  6  years.  Sensitivities  are  calculated  with 
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1  the  model  output  from  Aug.  1,  1966  to  Aug.  31,  1966.  This  process  is  looped  for  the 

2  two  wave-affected  parameters.  Fig.  A2  shows  the  ensemble  spread  of  T  with  respect 

3  to  a  and  P  at  different  depths.  The  ensemble  spread  of  T  near  the  sea  surface  is 

4  more  than  0.09  with  respect  to  /?  and  less  than  0.02  with  respect  to  a  .  The 

5  sensitivity  of  T  is  obviously  larger  to  /?  than  to  a  for  the  whole  depth, 

6  especially  in  the  upper  30m.  Small  sensitivity  in  the  lower  layer  indicates  that  the 

7  noise  may  be  stronger  than  the  signal  during  the  parameter  estimation  when  the  lower 

8  layer  temperature  observations  are  assimilated  into  the  bias  simulation  model. 

9  The  sensitivities  with  respect  to  the  wave-affected  parameters  are  also 

10  investigated  through  calculating  the  gradients  of  the  cost  function  with  the  parameters, 

qj  (P)J 

1 1  namely  —  and  —  .  Table  A1  shows  the  dependence  of  the  sensitivity  on  the  initial 

da  dp 

12  values  of  the  parameter  a  and  p .  When  the  initial  parameter  values  (a,  P)  are  set 

13  exactly  to  the  “truth”  values  (200,  2),  both  sensitivities  are  very  close  to  zero.  In 

14  general,  the  sensitivity  is  several  orders  of  magnitude  greater  on  P  than  on  a .  It 

15  indicates  that  the  parameter  a  is  more  vulnerable  to  be  disturbed  by  the  noises 

16  arisen  from  the  observational  errors  and  the  biased  initial  state  fields  during  the 

17  parameter  estimation. 

18  Appendix  B.  Correctness  test  of  the  gradient 

19  The  code  of  the  adjoint  model  is  produced  directly  through  the  Adjoint  Model 

20  Compiler  (AMC)  developed  by  Giering  and  Kaminski  (1998)  (Of  course,  a 

21  hand-coding  correction  is  necessary  after  that),  which  means  that  it  is  not  essential  to 

22  produce  the  code  of  the  tangent  linear  model  explicitly.  Therefore,  only  the 

23  correctness  of  the  adjoint  model  is  tested  here. 
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According  to  the  Taylor  expression,  one  has 


2  lim  cp{s)  =  lim  J(x°  gV/(Xo))  J(Xq)  x\  (A  1 ) 

«-*>  «-*>  _  s  <  W(x0),  W(x0)  > 

3  where  xo  is  any  control  variable,  the  symbol  <  >  represents  the  inner  product.  Fig. 

4  A3  shows  the  correctness  test  of  the  gradient  of  the  cost  function  with  respect  to  a 

5  and  /?  using  the  Eq.  (Al).  With  respect  to  a  ,  <p{s)  converges  to  1  as  £ 

6  decreases  from  10"3  to  10  s,  and  decreases  from  1  to  0.38  as  £  decreases  from  10  s 

7  to  10"10,  which  indicates  the  dominance  of  the  computational  errors  in  <p(s) .  With 

8  respect  to  p,  cp(s)  converges  to  1  as  s  decreases  from  10"6.  Therefore,  the  adjoint 

9  coding  is  valid. 
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1 


Figure  captions 


2  Figure  1.  Daily  temperature  in  1966  at  (a)  0  m,  (b)  10  m,  (c)  20  m  and  (d)  30  m  at 

3  the  OWS  Papa  with  the  “truth”  model  simulation  (solid  curve)  and  the  biased 

4  simulation  (dashed  curve). 

5  Figure  2.  Flowchart  of  the  wave-affected  parameter  estimation  using  the  variational 

6  method. 

7  Figure  3.  Time  series  of  the  estimated  wave-effected  parameters  (a)  a ,  and  (b)  /3 

8  for  PE  from  Aug.l  to  Aug.  30,  1966  (solid  curve),  where  both  assimilation 

9  window  and  frequency  are  1  day  and  the  depth  of  the  assimilation  is  30  m.  Here,  the 

10  dashes  curves  show  the  “truth”  (a,  p)  values. 

1 1  Figure  4.  Dependence  of  (a)  the  cost  function  and  (b)  the  norm  of  the  gradient  on  the 

12  number  of  iterations  on  Aug,  2,  1966. 

13  Figure  5.  Temporal  variations  of  the  natural  logarithm  of  the  cost  function  at  OWS 

14  Papa  from  Aug.  1  to  Aug.  30,  1966.  Here,  the  red,  blue  and  black  curves  represent  the 

15  total,  background,  and  observation  terms  of  the  cost  function. 

16  Figure  6.  Time  series  of  the  estimated  wave-effected  parameters  (a)  a ,  and  (b)  [3 

17  for  different  initial  parameter  values  from  Aug.l  to  Aug.  30,  1966.  Here,  both 

18  assimilation  window  and  frequency  are  1  day  and  the  depth  of  the  assimilation  is  30  m. 

19  The  black,  blue,  green,  yellow,  red,  pink,  purple,  orange  and  gray  solid  curves  in  panel 

20  (a)  and  the  corresponding  dashed  curves  in  panel  (b)  show  values  of  (  a  ,  (3)  with 

21  the  initial  guess  values  of  (0,0),  (100,2),  (100,3),  (200,1),  (200,3),  (300,1),  (300,2), 

22  (300,3)  and  (400,4),  respectively. 
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1  Figure  7.  Same  as  Figure  3  but  for  PE_SST  with  only  the  SST  observations  being 

2  assimilated. 

3  Figure  8.  Time  series  of  the  estimated  wave-effected  parameters  (a)  a ,  and  (b)  /? 

4  for  different  assimilation  window  and  frequency  from  Aug.l  to  Aug.  30,  1966  with 

5  30  m  as  the  depth  of  the  assimilation.  Here,  the  black,  blue,  red,  pink  solid  curves  in 

6  panel  (a)  and  panel  (b)  show  the  assimilation  frequency  are  48,  72,  96,  and  120  hours. 

7  Figure  9.  Temporal  variations  of  the  natural  logarithm  of  the  cost  function  at  OWS 

8  Papa  from  Aug.  1  to  Aug.  30,  1966.  The  solid  and  dashed  curves  represent  the  PE  and 

9  PE_5d  with  black  dots  denoting  the  time  that  the  observational  data  are  assimilated. 

10  Figure  10.  (a)  Sea  surface  temperature  and  (b)  mixed  layer  depth  from  CTRL  (black 

11  dashed  curve)  and  RPE  (black  solid  curve),  and  observations  (red  solid  curve)  at 

12  OWS  Papa  from  Aug.  1  to  Aug.  30,  1966.  The  horizontal  axis  represents  the  day 

13  relative  to  Aug.  1,  1966. 

14  Figure  11.  Sea  surface  temperature  from  observations  (red  solid  curve),  RSEPo 

15  (black  solid  curve),  and  RSE_Pd  (black  dashed  curve)  at  OWS  Papa  from  Aug.  1  to 

16  Sept  30,  1966.  The  horizontal  axis  represents  the  day  relative  to  Aug.  1,  1966. 

17  Figure  12.  Vertical  profiles  of  the  simulated  monthly-averaged  (a)  two  times  turbulent 

18  kinetic  energy  cf  (mV2),  and  (b)  and  vertical  mixing  coefficient  for  temperature  KH 

19  (10"3m2s_1)  from  RSE  Po  (solid  curve)  and  RSE  Pd  (dashed  curve)  at  OWS  Papa  in 

20  Sep.,  1966. 

21  Figure  Al.  Dependence  of  the  cost  function  on  a  and  p  for  (a)  10  >  P>  0,  and  (b)  3 

22  >  p  >  0. 
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1  Figure  A2.  Depth-dependence  of  ensemble  spread  of  temperature  with  respect  to  the 


2  wave-effected  parameters  a  (dashed  curve)  and  (3  (solid  curve). 

3  Figure  A3.  The  correctness  test  of  the  gradient  with  respect  to  (a)  a  ,  and  (b)  / 3 . 
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7 

8 


41 
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2  Table  1 .  All  assimilation  experiments  and  simulations  within  the  identical  synthetic 

3  experiment  framework. 


Name 

Description 

Control 

variables 

Assimilation 

windows 

Assimilation 

frequency 

Assimilation 

depth 

“Truth” 

model 

simulation 

O' =  200 

P  =  2 

- 

- 

- 

- 

Biased 

simulation 

a  =  100 

P=\ 

- 

- 

- 

- 

PE 

Parameter 

estimation 

rpn  rriYl—X 

? 

a  ,/3 

1  day 

1  day 

30  m 

PESST 

Parameter 

estimation 

rrin  rjm—\ 

a,P 

1  day 

1  day 

Sea  surface 

PE  [5  TT 

Parameter 

estimation  with 

the  “perfect” 

initial  fields 

derived  from 

the  “truth” 

P 

1  day 

1  day 

30  m 
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model 

simulation 

PE_/?_BI 

Parameter 

estimation  with 

the  “biased” 

initial  fields 

derived  from 

the  biased 

simulation 

P 

1  day 

1  day 

30  m 
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1  Table  2.  dependence  of  the  optimally-estimated  (a  ,J3)  on  the  standard  deviation  of 


2  temperature  observation 


Standard 

deviation  of 

temperature 

observation 

Estimated 

value  of  a 

Estimated 

value  of  /? 

Relative  error  of 

a 

Relative  error  of 

P 

icr3 

206.125 

1.982 

96.9% 

99.1% 

icr2 

167.343 

2.098 

83.6% 

95.1% 

0.05 

120.033 

2.114 

60.0% 

94.3% 

0.1 

100.096 

1.889 

50.0% 

94.5% 

0.5 

100.068 

0.866 

50.0% 

43.3% 
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1  Table  3.  Evolution  of  the  cost  function,  norm  of  the  projected  gradient  and  value  of 


2  ) 3  with  respect  to  the  number  of  iterations  for  the  direct  perturbed  method  with  the 

3  perfect  initial  field. 


Iteration  step 

Cost  function 

Norm  of  the 

projected 

gradient 

Value  of  f3 

0 

5.881 

2.097 

1.0 

1 

I.354e-5 

7.406e-2 

2.000365 

2 

2.63  le-9 

1.032e-3 

2.000005 

3 

7.44  le- 17 

3.042e-7 

2.000000 

4 

5.056e-17 

5.693e-9 

1.999999 

4 

5 
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Table  4.  Same  as  Table  3  but  with  the  biased  initial  field. 


Iteration  step 

Cost  function 

Norm  of  the 

projected 

gradient 

Value  of  /? 

0 

2.319e2 

4.819 

1.0 

1 

1.072e2 

3.351 

3.350811 

2 

1.071e2 

1.234 

3.317216 

3 

1.071e2 

9.003e-2 

3.301275 

4 

1.071e2 

1.982e-3 

3.302359 

5 

1.071e2 

6.043e-6 

3.302335 

6 

1.071e2 

3.627e-6 

3.302335 
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Table  5.  All  assimilation  experiments  and  model  simulations  within  the  real 


2  experiment  framework. 


Name 

Description 

Control 

variables 

Assimilation 

windows 

Assimilation 

frequency 

Assimilation 

depth 

Initial 

fields 

CTRL 

Simulation 

with 

a  =  200 

P =2 

- 

- 

- 

- 

1  Aug.  1966 

from  the  “truth” 

model  simulation 

RPE 

Real 

parameter 

estimation 

rj^n  rj-in—l 

a  ,p 

3  day 

3  day 

30  m 

Same  as  CTRL 

RSEPo 

Simulation 

using  the 

parameters 

estimated  by 

RPE 

- 

- 

- 

- 

Aug.  31,1966, 

derived  from 

RPE 

Simulation 

using  the 

same 

parameters 

as  in  CTRL 

- 

- 

- 

- 

Same  as  RSE  Po 

Table  6.  Evolution  of  the  cost  function,  (a  ,  ft)  with  respect  to  the  number  of 
iterations  for  the  real  assimilation. 

Iteration  step  Normalized  Value  of  a  Value  of  f3 
cost  function 


1  Table  A1  Dependence  of  the  sensitivity  on  the  initial  values  of  the  parameter  a  and  j3 


2 


3 

4 


Initial  values  of 

(a,j3) 

Sensitivity  of  a 

Sensitivity  of  J3 

(0,0) 

-7.5x1 0"5 

2.6x1 04 

(100,1) 

-4.69 

-574.96 

(100,2) 

-42.63 

-5325.69 

(100,3) 

158.70 

1.79x104 

(200,1) 

-7.41 

-4427.36 

(200,2) 

4.0x10-11 

-3.3x10-10 

(200,3) 

172.36 

2.73x104 

(300,1) 

-23.82 

-1.61x104 

(300,2) 

7.76 

2.50x103 

(300,3) 

429.53 

1.57x105 

(400,4) 

471.18 

1.80x10s 
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1 


3  Figurel.  Daily  temperature  in  1966  at  (a)  Om,  (b)  10m,  (c)  20m  and  (d)  30m  at  the 

4  OWS  Papa  with  the  “truth”  model  simulation  (solid  curve)  and  the  biased  simulation 

5  (dashed  curve). 
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1 

2 


26  Figure  2.  Flowchart  of  the  wave-affected  parameter  estimation  with  the  variational 

27  method. 

28 
29 
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Beta  A|Pha 


200 


3  for  PE  from  Aug.l  to  Aug.  30,  1966  (solid  curve),  where  both  assimilation  window 

4  and  frequency  are  1  day  and  the  depth  of  the  assimilation  is  30m.  The  dashes  curves 

5  show  the  “truth”  (a,  p)  values. 

6 
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3  number  of  iterations  on  Aug,  2,  1966. 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 


Papa  from  Aug.  1  to  Aug.  30,  1966.  Red,  blue  and  black  lines  are  the  total  terms,  the 
background  term  and  the  observation  term  of  the  cost  function  respectively. 
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2  Figure  6.  Time  series  of  the  estimated  wave-effected  parameters  (a)  a ,  and  (b)  /? 

3  for  different  initial  parameter  values  from  Aug.l  to  Aug.  30,  1966,  where  both 

4  assimilation  window  and  frequency  are  1  day  and  the  depth  of  the  assimilation  is  30m. 

5  Black,  blue,  green,  yellow,  red,  pink,  purple,  orange  and  gray  solid  lines  in  panel  (a) 

6  and  the  corresponding  dash  lines  in  panel  (b)  show  values  of  ( a  ,  /?)  =  (0,0),  (100,2), 

7  (100,3),  (200,1),  (200,3),  (300,1),  (300,2),  (300,3)  and  (400,4)  respectively. 
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2  Figure  7.  Same  as  Figure  3,  but  for  PE_SST,  where  only  the  SST  observations  are 

3  assimilated. 
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240 


3  for  different  assimilation  window  and  frequency  from  Aug.l  to  Aug.  30,  1966,  where 

4  the  depth  of  the  assimilation  is  30m.  Black,  blue,  red,  pink  solid  lines  in  panel  (a)  and 

5  panel  (b)  show  the  assimilation  frequency  are  48,  72,  96  and  120  hours. 

6  . 


7 


57 


3  Papa  from  Aug.  1  to  Aug.  30,  1966.  The  solid  and  dashed  curves  represent  PE  and 

4  PE_5d  with  black  dots  denoting  the  time  that  observations  are  assimilated. 
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3  Figure  10.  (a)  Sea  surface  temperature  and  (b)  mixed  layer  depth  from  CTRL  (black 

4  dashed  curve)  and  RPE  (black  solid  curve),  and  observations  (red  solid  curve)  at 

5  OWS  Papa  from  Aug.  1  to  Aug.  30,  1966.  The  horizontal  axis  represents  the  day 

6  relative  to  Aug.  1,  1966. 
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3  (black  solid  curve),  and  RSE  Pd  (black  dashed  curve)  at  OWS  Papa  from  Aug.  1  to 

4  Sept  30,  1966.  The  horizontal  axis  represents  the  day  relative  to  Aug.  1,  1966. 
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2  Figure  12.  Vertical  profiles  of  the  simulated  monthly-averaged  (a)  two  times  turbulent 

3  kinetic  energy  (mV2 3 4 5 6),  and  (b)  and  vertical  mixing  coefficient  for  temperature 

4  (KTmV1)  from  RSEPo  (solid  curve)  and  RSE_Pd  (dashed  curve)  at  OWS  Papa  in 

5  Sept.,  1966. 
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Figure  Al.  Dependence  of  the  cost  function  on  a  and  P  with  (a)  10  >  P>  0,  and  (b)  3  > 


P>0. 
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3  parameters  a  (dashed  curve)  and  /?  (solid  curve)  at  different  depths. 
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Figure  A3.  The  correctness  test  of  the  gradient  with  respect  to  (a)  a  ,  and  (b)  f:S . 
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